OISST Sea Surface Temperature Report
Adam A. Kemberling
Updated on: 2021-02-17
Gulf of Maine Sea Surface Temperature Outlook
Region Extent
The spatial extent for Gulf of Maine is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.
# File paths for various extents based on params$region
poly_path <- switch(
EXPR = tolower(params$region),
"gulf of maine" = "Shapefiles/gmri_sst_focal_areas/apershing_gulf_of_maine.geojson",
"cpr gulf of maine" = "Shapefiles/gmri_sst_focal_areas/cpr_gulf_of_maine.geojson",
"northwest atlantic" = "Shapefiles/gmri_sst_focal_areas/aak_northwest_atlantic.geojson")
# Load the bounding box for Andy's GOM to show they align
poly_path <- paste0(res_path, poly_path)
region_extent <- st_read(poly_path, quiet = TRUE)
# Pull extents for the region for crop
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]
# Zoom out for cpr extent
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-70.875, -65.375)
crop_y <- c(40.375, 45.125)}
# Full plot
ggplot() +
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
geom_sf(data = region_extent, color = "royalblue",
fill = "transparent", linetype = 2, size = 1) +
map_theme +
coord_sf(xlim = crop_x,
ylim = crop_y, expand = T) 
GMRI OISST Resources on Box
For many of our most-used areas, regional climatologies and anomalies have been processed using a handful of jupyter notebook work flows to take advantage of {xarray}. The following sections will highlight what those work flows are, what they detail, and how they can be updated using the Gulf of Maine as an example.
The root folder for the following OISST products is ~/Box/RES_Data/OISST/oisst_mainstays.
Regional Timeseries
A region-masked time series is the most basic building block for relaying information for any of these areas. For any desired area (represented by a spatial polygon) a time series table of observed sea surface temperature is returned that contains the mean sea surface temperature for that area only for each time step.
Rather than repeat this process again for the climatology table and the anomalies these three values and a few additional fields are stored in the same table.
Pre-processed tables have been placed on box as part of the NSF C-ACCEL project under the folder: ~/Box/RES_Data/OISST/oisst_mainstays/regional_timeseries. A function to quickly access these timelines has been added to the {gmRi} package as oisst_access_timeseries().
# Switch path based on local vs. docker
timeseries_path <- switch(
EXPR = tolower(params$region),
"gulf of maine" = "OISSTv2_anom_apershing_gulf_of_maine.csv",
"cpr gulf of maine" = "OISSTv2_anom_cpr_gulf_of_maine.csv",
"northwest atlantic" = "OISSTv2_anom_northwest_atlantic.csv")
# Append to form complete path
timeseries_path <- str_c(
oisst_path,
"regional_timeseries/gmri_sst_focal_areas/",
timeseries_path)
# Load timeseries
region_timeseries <- read_csv(timeseries_path,
guess_max = 1e5,
col_types = cols())
# Display Table of first 6 entries
head(region_timeseries) %>%
mutate_if(is_numeric, round, 2) %>%
kable() %>%
kable_styling()| time | modified_ordinal_day | sst | sst_clim | sst_anom | clim_sd | log_lik |
|---|---|---|---|---|---|---|
| 1981-09-01 | 245 | 15.78 | 16.72 | -0.94 | 2.32 | 1.84 |
| 1981-09-02 | 246 | 15.79 | 16.67 | -0.89 | 2.29 | 1.82 |
| 1981-09-03 | 247 | 15.49 | 16.64 | -1.14 | 2.28 | 1.87 |
| 1981-09-04 | 248 | 14.99 | 16.60 | -1.60 | 2.25 | 1.98 |
| 1981-09-05 | 249 | 14.84 | 16.51 | -1.67 | 2.21 | 2.00 |
| 1981-09-06 | 250 | 15.08 | 16.49 | -1.41 | 2.18 | 1.91 |
Climatologies are currently set up to calculate daily averages on a modified julian day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not. This preserves comparisons across calendar dates such-as: “The average temperature on march 1st is ____ for the reference period ____ to ____”
In these tables sst is the mean temperature observed for that date averaged across all cells within the area. sst_clim & clim_sd are the climate means and standard deviations for a 1982-2011 climatology. sst_anom is the daily observed minus the climate mean, and log_lik is the log likelihood of getting that value given a normal distribution with mean of sst_clim and standard deviation of clim_sd.
Marine Heatwaves
The {heatwaveR} package provides a relatively quick way of working with tabular data to calculate a seasonal climate mean, as well as heatwaves and coldwaves at a desired threshold.
These functions can be used to get the climatology using the standard day of year if desired, and track the number and length of heatwave events. Heatwave events follow the definition of Hobday et al. 2016 and are set up to use a 90% threshold for heatwave/coldwave detection.
The heatwave history for the Gulf of Maine trawl region displayed above is as follows:
# Wrapper function to do heatwaves and coldwaves simultaneously at 90%
pull_heatwave_events <- function(temperature_timeseries) {
# Pull the two column dataframe for mhw estimation
test_ts <- data.frame(t = temperature_timeseries$time,
temp = temperature_timeseries$sst)
# Detect the events in a time series
ts <- ts2clm(data = test_ts, climatologyPeriod = c("1982-01-01", "2011-12-31"))
#heatwaves
mhw <- detect_event(ts)
# prep heatwave data
mhw_out <- mhw$climatology %>%
mutate(sst_anom = temp - seas) %>%
select(time = t,
sst = temp,
seas,
sst_anom,
mhw_thresh = thresh,
mhw_event = event,
mhw_event_no = event_no)
# 2. Detect cold spells
ts <- ts2clm(data = test_ts,
climatologyPeriod = c("1982-01-01", "2011-12-31"),
pctile = 10)
mcs <- detect_event(ts, coldSpells = TRUE) #coldSpells = TRUE flips boolean to < thresh
# prep cold spell data
mcs_out <- mcs$climatology %>%
select(time = t,
mcs_thresh = thresh,
mcs_event = event,
mcs_event_no = event_no)
# join heatwaves to coldwaves
hot_and_cold <- left_join(mhw_out, mcs_out, by = "time")
# 3. Data formatting for plotting,
# adds columns to plot hw and cs seperately
events_out <- hot_and_cold %>%
mutate(status = ifelse(mhw_event == TRUE, "Marine Heatwave Event", "Sea Surface Temperature"),
status = ifelse(mcs_event == TRUE, "Marine Cold Spell Event", status),
hwe = ifelse(mhw_event == TRUE, sst, NA),
cse = ifelse(mcs_event == TRUE, sst, NA),
nonevent = ifelse(mhw_event == FALSE & mcs_event == FALSE, sst, NA))
# Close the gaps between a mhw event and sst (might not need if full line for temp exists)
events_out <- events_out %>%
mutate(hwe = ifelse(is.na(hwe) & is.na(lag(hwe)) == FALSE, sst, hwe),
cse = ifelse(is.na(cse) & is.na(lag(cse)) == FALSE, sst, cse))
return(events_out)
}
# Use function to process heatwave data for plotting
region_heatwaves <- pull_heatwave_events(region_timeseries)Heatwave timelines can be then be plotted using the {ggplot2} package for static plots, or by using the {plotly} package for interactive content.
Interactive SST Timeline
For anything we wish to host on the web there is an option to display tables and graphs that are interactive. The {plotly} package is one-such tool for producing plots that allow users to pan, zoom, and highlight discrete observations.
# Helper function
plotly_mhw_plots <- function(data){
# How to get rgb() from colors
plot_cols <- list(
"gray20" = col2rgb(col = "gray20"),
"gray40" = col2rgb(col = "gray40"),
"royalblue" = col2rgb(col = "royalblue"),
"darkred" = col2rgb(col = "darkred"),
"royalblue" = col2rgb(col = "royalblue"))
# Building the plot
fig <- plot_ly(data, x = ~time)
# Sea Surface Temperature
fig <- fig %>% add_trace(y = ~sst,
name = 'Sea Surface Temperature',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(65, 105, 225)",
width = 2))
# Heatwave Threshold
fig <- fig %>% add_trace(y = ~mhw_thresh,
name = 'MHW Threshold ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(51, 51, 51)",
width = 1,
dash = 'dot'))
# Seasonal Climatology
fig <- fig %>% add_trace(y = ~seas,
name = 'Seasonal Climatology ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(102, 102, 102)",
width = 3,
dash = 'dash'))
# Marine Cold Spell Threshold
fig <- fig %>% add_trace(y = ~mcs_thresh,
name = 'MCS Threshold ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(51, 51, 51)",
width = 1,
dash = 'dot'))
# Heatwave Event
fig <- fig %>% add_trace(y = ~hwe,
name = 'Marine Heatwave Event ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(139, 0, 0)",
width = 2))
# Cold Spell Event
fig <- fig %>% add_trace(y = ~cse,
name = 'Marine Cold Spell Event ',
mode = 'lines',
type = "scatter",
line = list(color = "rgb(65, 105, 225)",
width = 2))
# Axis Formatting
fig <- fig %>% layout(xaxis = list(title = ""),
yaxis = list (title = "Temperature (degrees C)"))
# Legend formatting
fig <- fig %>% layout(legend = list(orientation = 'h'))
return(fig)
}
# Grab data from the most recent year through present day
last_year <- Sys.Date() - 365
last_year <- last_year - yday(last_year) + 1
last_yr_heatwaves <- region_heatwaves %>%
filter(time >= last_year)
# Plot that timeseries
last_yr_heatwaves %>%
filter(time >= last_year) %>%
plotly_mhw_plots()Static SST Timeline
In most cases a static image will suffice, in which case {ggplot2} is a familiar go-to.
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"Cold Spell Event" = "lightblue",
"MHW Threshold" = "gray30",
"MCS Threshold" = "gray30",
"Seasonal Climatology" = "gray30")
# Set the label with degree symbol
ylab <- expression("Sea Surface Temperature"~degree~C)
# How many heatwave events:
num_events <- max(last_yr_heatwaves$mhw_event_no, na.rm = T) - min(last_yr_heatwaves$mhw_event_no, na.rm = T)
# How many heatwave days
num_hw_days <- sum(last_yr_heatwaves$mhw_event, na.rm = T)
# Plot the last 365 days
ggplot(last_yr_heatwaves, aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = cse, color = "Cold Spell Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 2, size = .5) +
geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 2, size = .5) +
geom_line(aes(y = seas, color = "Seasonal Climatology"), size = .5) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
theme(legend.title = element_blank(),
legend.position = "bottom") +
labs(x = "", y = ylab,
caption = paste0("X-axis start date: ", last_year,
"\nNumber of Heatwave events: ", num_events,
"\nNumber of Heatwave Days: ", num_hw_days))
Warming Rates
Regional warming trends below were calculated using annual averages, and consequently the model parameters reflect annual increases in sea surface temperature.
Things to add:
* Global Trends (grey)
* Regional Trend (blue) * Coefficients on Annual scale
Annual
# Summarise by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>%
mutate(year = year(time)) %>%
group_by(year) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
.groups = "keep") %>% ungroup() %>%
mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))
# Plot the warming rate for the region
ggplot(annual_summary, aes(yr_as_dtime, sst_anom)) +
# Add daily data
geom_line(data = region_timeseries, aes(time, sst_anom),
alpha = 0.5, color = "gray") +
# Overlay yearly means
geom_line(alpha = 0.7) +
geom_point(alpha = 0.7) +
# Add regression line - all data
geom_smooth(method = "lm",
aes(color = "1982-2020 Regional Trend"),
formula = y ~ x, se = F,
linetype = 2) +
# Add equation
stat_poly_eq(formula = y ~ x,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~~")),
color = gmri_cols("gmri blue"),
parse = T) +
# Add regression line - all data
geom_smooth(data = filter(annual_summary, year %in% c(2006:2020)),
method = "lm",
aes(color = "2006-2020 Regional Trend"),
formula = y ~ x, se = F,
linetype = 2) +
# Add equation
stat_poly_eq(data = filter(annual_summary, year %in% c(2006:2020)),
formula = y ~ x,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~~")),
color = gmri_cols("orange"),
parse = T, label.x = 0.05, label.y = 0.9) +
# Colors, labels, theme
scale_color_manual(values = c(
"1982-2020 Regional Trend" = as.character(gmri_cols("gmri blue")),
"2006-2020 Regional Trend" = as.character(gmri_cols("orange")))) +
labs(x = "", y = "Sea Surface Temperature Anomaly",
caption = "Regression Coefficients Reflect Annual Change in Degree C.") +
theme(legend.title = element_blank(),
legend.position = c(0.85, 0.1),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent"))
Quarterly
# Repeat for the seasons (equal quarters here*)
quarter_summary <- region_timeseries %>%
mutate(year = year(time),
season = factor(quarter(time, fiscal_start = 1)),
season = fct_recode(season, c("Q1" = "1"), c("Q2" = "2"),
c("Q3" = "3"), c("Q4" = "4"))) %>%
group_by(year, season) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
.groups = "keep")
quarter_summary %>%
ggplot(aes(year, sst_anom)) +
geom_line(group = 1) +
geom_point() +
geom_smooth(method = "lm",
aes(color = "Regional Trend"),
formula = y ~ x, se = F, linetype = 2) +
stat_poly_eq(formula = y ~ x,
color = gmri_cols("gmri blue"),
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = T) +
scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
labs(x = "", y = "Sea Surface Temperature Anomaly") +
theme(legend.title = element_blank(),
legend.position = c(0.875, 0.05),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent")) +
facet_wrap(~season)
Heatwave Events
As mentioned above heatwave events were determined using the methods of Hobday et al. 2016 and implemented using {heatwaveR} with a 90% threshold. The region used is again that same polygon outlined in the first figure, the bottom right panel of the maps at the beginning of the document.
# Prep the legend title
guide_lab <- expression("Sea Surface Temperature Anomaly"~degree~C)
# get new axis dimensions, y = year, x = day within year
# use flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_heatwaves %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# assemble plot
ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2020.5), expand = c(0,0)) +
labs(x = "", y = "") +
#scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme_classic() +
theme(legend.title = element_text(angle = 90))
Anomaly Maps
For demonstration purposes I will only be loading the 2020 global anomaly data for the below images.
# Access information to netcdf on box
nc_year <- "2020"
anom_path <- str_c(oisst_path, "annual_anomalies/daily_anoms_", nc_year, ".nc")
# Load 2020 as stack
anoms_2020 <- stack(anom_path)Mean Global
The following code will subset the anomalies for July and plot the average sea surface temperature anomalies for that month:
# Get the mean temperature anomalies for July
july_dates <- which(str_sub(names(anoms_2020), 7, 8) == "07")
july_avg <- mean(anoms_2020[[july_dates]])
# Convert wgs84 raster to stars array
july_st <- st_as_stars(rotate(july_avg))
# Plot global map
ggplot() +
geom_stars(data = july_st) +
geom_sf(data = world, fill = "gray90") +
# scale_fill_gradient2(low = "blue",
# mid = "white",
# high = "red") +
scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
map_theme +
coord_sf(expand = FALSE) +
guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
Regional Map
Using the July average sst anomaly we can then clip down to the Gulf of Maine for an aerial view of the region. For this execution the bounding box of lat-lon coordinates from the regional extent will be used.
# Clip Raster - Convert to stars
shape_extent <- c(crop_x, crop_y)
region_ras <- crop(rotate(july_avg), extent(shape_extent))
region_st <- st_as_stars(region_ras)
# Get crop bounds for coord_sf
crop_x <- st_bbox(region_st)[c(1,3)]
crop_y <- st_bbox(region_st)[c(2,4)]
# Zoom out for cpr extent, same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-71, -65.5)
crop_y <- c(40.5, 45)}
# Plot Gulf of Maine - wgs84
ggplot() +
geom_stars(data = region_st) +
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
map_theme +
coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
guides("fill" = guide_colorbar(
title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
Regionally Projected Map
These lat/lon region extent boxes are good showcases of how something that is “rectangular” by its lat/lon dimensions that has its shape change to when projected onto a non-flat surface of the earth.
To get the OISST anomaly grid to display on a curviliniear projection the originals were cropped using the corner coordinates of the extent plotted above. The raster was then warped into a regular grid using an appropriate projection as displayed below.
For the most dramatic case set the region parameter to “Northwest Atlantic”
#### Candidate coordinate reference systems ####
#### NOTE:
# transformation from rasters different than polygons
# as cells need to stretch and bend ####
# Resource: Tranformation vs. Warping
# https://r-spatial.github.io/stars/articles/stars5.html
# Robinson projection
robinson_proj <- "+proj=robin"
# Albers equal area: centered on -70 degrees
# The settings for lat_1 and lat_2 are the locations at which
# the cone intersects the earth, so distortion is minimized at those latitudes
alb_70 <- "+proj=aea +lat_1=30 +lat_2=50 +lon_0=-70"
# custom sterographic projection for this nw atlantic, centered using lon_0
stereographic_north <- "+proj=stere +lat_0=90 +lat_ts=75 +lon_0=-57"
# equal earth projection - not working, don't think sf has functionality for it yet
eqearth_proj <- "+proj=eqearth"
# Choose crs using parameter
projection_crs <- switch(
tolower(params$region),
"gulf of maine" = alb_70,
"cpr gulf of maine" = alb_70,
"northwest atlantic" = stereographic_north)
# Transform all the polygons
region_projected <- st_transform(region_extent, crs = projection_crs)
canada_projected <- st_transform(canada, crs = projection_crs)
newengland_projected <- st_transform(new_england, crs = projection_crs)
greenland_projected <- st_transform(greenland, crs = projection_crs)
# coord_sf Crop bounds in projection units for coord_sf
crop_x <- st_bbox(region_projected)[c(1,3)]
crop_y <- st_bbox(region_projected)[c(2,4)]
# Zoom out for cpr extent, same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-73185.13, 386673.34)
crop_y <- c(4269044, 4813146)}
# Lower the ymin a touch for NW Atlantic
if(tolower(params$region) == "northwest atlantic"){crop_y <- crop_y - c(100000, 0)}
# Warp to grid of chosen CRS
projection_grid <- region_st %>%
st_transform(projection_crs) %>%
st_bbox() %>%
st_as_stars()
region_warp_ras <- region_st %>%
st_warp(projection_grid)
# Plot everything together
ggplot() +
geom_stars(data = region_warp_ras) +
geom_sf(data = newengland_projected, fill = "gray90") +
geom_sf(data = canada_projected, fill = "gray90") +
geom_sf(data = greenland_projected, fill = "gray90") +
coord_sf(crs = projection_crs,
xlim = crop_x, ylim = crop_y, expand = T) +
map_theme +
scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
theme(
panel.border = element_rect(color = "black", fill = NA),
plot.background = element_rect(color = "transparent", fill = "transparent"),
axis.title.x = element_blank(), # turn off titles
axis.title.y = element_blank(),
legend.position = "bottom",
legend.title.align = 0.5)
Warming Rates
Global warming rates are calculated by getting the average temperature for each year across all cells in the data, then using a linear regression for each to determine the annual warming rate in degrees Celsius.
The warming rates are then ranked from high to low to identify areas with the most rapid ocean warming. The following maps display the warming rates and the percentile ranks for the top 20% of warming rates, or rates above the 80th percentile of the data.
# Access information to warming rate netcdf files on box
rates_path <- paste0(oisst_path, "warming_rates/annual_warming_rates")
# Going to reclassify raster from continuous to discrete:
# Source: https://www.earthdatascience.org/courses/earth-analytics/lidar-raster-data-r/classify-raster/
reclassify_to_discrete <- function(ranks_stack, rates_stack, percentile_cutoff = 80){
# create classification matrix
reclass_df <- c(0.00, 0.05, 0, #####
0.05, 0.10, 5,
0.10, 0.15, 10,
0.15, 0.20, 15,
0.20, 0.25, 20,
0.25, 0.30, 25,
0.30, 0.35, 30,
0.35, 0.40, 35,
0.40, 0.45, 40,
0.45, 0.50, 45,
0.50, 0.55, 50,
0.55, 0.60, 55,
0.60, 0.65, 60,
0.65, 0.70, 65,
0.70, 0.75, 70,
0.75, 0.80, 75,
0.80, 0.85, 80,
0.85, 0.90, 85,
0.90, 0.95, 90,
0.95, 1, 95)
#####
# reshape the object into a matrix with columns and rows
reclass_m <- matrix(reclass_df,
ncol = 3,
byrow = TRUE)
# reclassify the raster using the reclass object - reclass_m
ranks_classified <- reclassify(ranks_stack,
reclass_m)
# Save the un-masked layers as stars objects
rates_raw_st <- st_as_stars(rotate(rates_stack))
ranks_raw_st <- st_as_stars(rotate(ranks_stack))
# Masking Below percentile cutoff - for ranking raster and rates raster
ranks_stack[ranks_classified < percentile_cutoff] <- NA
rates_stack[ranks_classified < percentile_cutoff] <- NA
# Converting to stars
rates_st <- st_as_stars(rotate(rates_stack))
ranks_st <- st_as_stars(rotate(ranks_stack))
ranks_c_st <- st_as_stars(rotate(ranks_classified))
# #get scale ranges so they still pop
# rate_range <- c("min" = cellStats(rates_stack, 'min')[1],
# "max" = cellStats(rates_stack,'max')[1])
# Return the three stars objects
return(list("rates_raw" = rates_raw_st,
"ranks_raw" = ranks_raw_st,
"rates" = rates_st,
"ranks" = ranks_st,
"ranks_discrete" = ranks_c_st))
}1982-2011
#### Data Prep ####
percentile_cutoff <- 80
# 1982-2020
rates_stack_all <- stack(str_c(rates_path, "1982to2020.nc"),
varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2020.nc"),
varname = "rate_percentile")
rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all,
ranks_stack = ranks_stack_all,
percentile_cutoff = percentile_cutoff)
rates_st <- rates_prepped$rates
rates_raw_st <- rates_prepped$rates_raw
ranks_st <- rates_prepped$ranksGlobal Rates
#### Plot Setup ####
# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)
# ggplot using stars
ggplot() +
geom_stars(data = rates_st) +
#geom_stars(data = rates_raw_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 1982-2020"))
Global Ranks
# ggplot using stars
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 1982-2020"))
2006-2011
#### Data Prep ####
# 2006-2020
rates_stack_all <- stack(str_c(rates_path, "2006to2020.nc"),
varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "2006to2020.nc"),
varname = "rate_percentile")
rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all,
ranks_stack = ranks_stack_all)
rates_st <- rates_prepped$rates
ranks_st <- rates_prepped$ranksGlobal Rates
#### Plot Setup ####
# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)
# ggplot using stars
ggplot() +
geom_stars(data = rates_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 2006-2020"))
Global Ranks
# ggplot using stars
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 2006-2020"))
A work by Adam A. Kemberling
Akemberling@gmri.org